4 Psykiska/psykosomatiska besvär

4.1 Bakgrund

Item/frågor har etiketter F88-F99 i datafilen, och motsvaras av fråga 90-101 i PDF-filen med frågor.

Samtliga frågor har fem svarskategorier, vilka varierar mellan frågorna. Fem frågor har svarskategorier från “Aldrig” till “Flera gånger i veckan”. Sex frågor har från “Sällan” till “Väldigt ofta”, och en från “Inte alls” till “Väldigt mycket”.

Svarsdata har kodats så att högre poäng innebär mera besvär/högre risk.

Sektionen i enkäten inleds med meningen: “NÅGRA FRÅGOR OM HUR DU MÅR”.

4.2 Lista med enkätfrågorna

Code

RIlistitems(df.omit.na)

itemnr	item
F88	Hur ofta har du haft huvudvärk detta läsår?
F89	Känner du dig ledsen och deppig utan att veta varför?
F90	Händer det att du känner dig rädd utan att veta varför?
F91	Hur ofta har du dålig aptit?
F92	Hur mycket skulle du vilja ändra på dig själv?
F93	Hur ofta har du under detta läsår haft ”nervös mage” (t.ex. magknip, magkramper, orolig mage, illamående, gaser, förstoppning eller diarré)?
F94	Hur ofta tycker du att du inget duger till?
F95	Hur ofta har du under detta läsår haft svårt att somna?
F96	Är du nöjd med ditt utseende?
F97	Känner du dig slö och olustig?
F98	Hur ofta har det hänt under detta läsår att du sovit oroligt och vaknat under natten?
F99	Hur ofta tycker du att det är riktigt härligt att leva?

4.2.1 Demografi

Vi har 9595 deltagare i samplet från 2014, och deras könsfördelning återges i tabellen nedan. Deltagare som saknar data på samtliga frågor är borttagna ur analysen.

4.3 Deskriptiva data

4.3.1 Demografi

Code

RIdemographics(dif.gender, "Kön")
RIdemographics(dif.arskurs, "Årskurs")

Kön	n	Percent
Flicka	5179	54
Pojke	4416	46

Årskurs	n	Percent
Åk 9	4109	42.8
Gy 2	5486	57.2

4.3.2 Item-data

Tile plot

Code

RItileplot(df.omit.na)

Stacked bars

Code

RIbarstack(df.omit.na)

Barplots

Code

RIbarplot(df.omit.na)

4.4 Rasch-analys 1 samtliga items

itemnr	item
F88	Hur ofta har du haft huvudvärk detta läsår?
F89	Känner du dig ledsen och deppig utan att veta varför?
F90	Händer det att du känner dig rädd utan att veta varför?
F91	Hur ofta har du dålig aptit?
F92	Hur mycket skulle du vilja ändra på dig själv?
F93	Hur ofta har du under detta läsår haft ”nervös mage” (t.ex. magknip, magkramper, orolig mage, illamående, gaser, förstoppning eller diarré)?
F94	Hur ofta tycker du att du inget duger till?
F95	Hur ofta har du under detta läsår haft svårt att somna?
F96	Är du nöjd med ditt utseende?
F97	Känner du dig slö och olustig?
F98	Hur ofta har det hänt under detta läsår att du sovit oroligt och vaknat under natten?
F99	Hur ofta tycker du att det är riktigt härligt att leva?

Item fit

Code

RIitemfitPCM2(df.omit.na, 300, 32, 8)

	OutfitMSQ	InfitMSQ	OutfitZSTD	InfitZSTD
F88	1.06	1.043	0.749	0.775
F89	0.709	0.716	-4.108	-4.521
F90	0.851	0.896	-0.773	-1.232
F91	1.133	1.094	1.129	1.354
F92	0.788	0.794	-2.92	-3.126
F93	1.029	1.015	0.405	0.38
F94	0.839	0.819	-2.234	-2.329
F95	1.059	1.054	0.924	0.694
F96	0.995	0.969	0.27	-0.508
F97	0.968	0.959	-0.389	-0.415
F98	0.956	0.965	-0.557	-0.561
F99	1.072	0.995	0.507	-0.099

PCA

Code

RIpcmPCA(df.omit.na)

PCA of Rasch model residuals

Eigenvalues
2.02
1.31
1.22
1.11
1.08

Loadings 1st contrast

Code

RIloadLoc(df.omit.na)

Residualkorrelationer

Code

RIresidcorr(df.omit.na, cutoff = 0.2)

	F88	F89	F90	F91	F92	F93	F94	F95	F96	F97	F98	F99
F88
F89	-0.09
F90	-0.06	0.11
F91	-0.05	-0.13	-0.1
F92	-0.2	-0.01	-0.08	-0.12
F93	0.02	-0.08	-0.06	0	-0.1
F94	-0.18	0.06	-0.03	-0.21	0.12	-0.15
F95	-0.02	-0.17	-0.11	-0.02	-0.2	-0.06	-0.24
F96	-0.19	-0.12	-0.12	-0.22	0.29	-0.21	0.17	-0.24
F97	-0.11	-0.02	-0.05	-0.16	-0.07	-0.17	0.02	-0.13	-0.01
F98	0.02	-0.11	-0.04	-0.07	-0.21	-0.06	-0.17	0.18	-0.22	-0.14
F99	-0.19	0.02	-0.1	-0.13	-0.01	-0.17	-0.03	-0.13	0.08	-0.02	-0.17
Note:
Relative cut-off value (highlighted in red) is 0.122, which is 0.2 above the average correlation.

Targeting

Code

RItargeting(df.omit.na)

Itemhierarki

Code

RIitemHierarchy(df.omit.na)

Svarskategorier

Code

plot(mirt.rasch, type="trace")

PCA på residualer från Rasch-modellen indikerar möjlig multidimensionalitet. Utifrån faktorladdningarna på första residualkontrasten ser vi två potentiella kluster i data. Det ena betecknas av de psykosomatiska frågorna, det andra av psykiska besvär. Vi tittar på klustren separat.

4.5 Rasch-analys 1 - psykosomatiska besvär

Code

items.psm <- c("F88","F91","F93","F95","F98")

df.psm <- df.omit.na %>% 
  select(any_of(items.psm))

RIlistItemsMargin(df.psm, fontsize = 13)

itemnr	item
F88	Hur ofta har du haft huvudvärk detta läsår?
F91	Hur ofta har du dålig aptit?
F93	Hur ofta har du under detta läsår haft ”nervös mage” (t.ex. magknip, magkramper, orolig mage, illamående, gaser, förstoppning eller diarré)?
F95	Hur ofta har du under detta läsår haft svårt att somna?
F98	Hur ofta har det hänt under detta läsår att du sovit oroligt och vaknat under natten?

Item fit

Code

RIitemfitPCM2(df.psm, 300, 32, 8)

	OutfitMSQ	InfitMSQ	OutfitZSTD	InfitZSTD
F88	0.897	0.894	-1.126	-1.216
F91	0.908	0.919	-0.996	-0.98
F93	0.894	0.889	-1.127	-1.712
F95	0.783	0.802	-2.589	-2.966
F98	0.77	0.79	-3.037	-2.947

PCA

Code

RIpcmPCA(df.psm)

PCA of Rasch model residuals

Eigenvalues
1.50
1.28
1.17
1.04
0.01

Residualkorrelationer

Code

RIresidcorr(df.psm, cutoff = 0.2)

	F88	F91	F93	F95	F98
F88
F91	-0.23
F93	-0.13	-0.17
F95	-0.22	-0.24	-0.26
F98	-0.14	-0.28	-0.23	0.02
Note:
Relative cut-off value (highlighted in red) is 0.012, which is 0.2 above the average correlation.

Targeting

Code

RItargeting(df.psm)

Itemhierarki

Code

RIitemHierarchy(df.psm)

Svarskategorier

Code

plot(mirt.rasch, type="trace")

Barplots

Code

for (i in 1:ncol(df.psm)) {
    barplot(table(df.psm[, i]), col = "#8dc8c7", main = names(df.psm[i]), 
      ylab = "Number of responses", xlab = (itemlabels %>% 
                                              filter(itemnr %in% names(df.psm)) 
                                            %>% .[i,2]))
  }

Bortsett från svarskategorierna fungerar items acceptabelt.

4.5.1 Omkodning av svarskategorier

Vi behöver åtgärda flera items genom att slå samman svarskategorier:

• F91: 1+2 och 3+4
• F93 och F98: 3+4
• F95: 0+1

Code

df.psm$F91 <- recode(df.psm$F91, "2=1;3=2;4=2", as.factor = FALSE)
df.psm$F93 <- recode(df.psm$F93, "4=3", as.factor = FALSE)
df.psm$F98 <- recode(df.psm$F98, "4=3", as.factor = FALSE)
df.psm$F95 <- recode(df.psm$F95, "4=3;3=2;2=1;1=0", as.factor = FALSE)

Svarskategorier

Code

plot(mirt.rasch, type="trace")

Itemhierarki (med konfidensintervall runt svarströsklar)

Code

RIitemHierarchy(df.psm)

4.6 Rasch-analys 2 - psykosomatiska besvär

Efter att ha åtgärdat svarskategorierna.

itemnr	item
F88	Hur ofta har du haft huvudvärk detta läsår?
F91	Hur ofta har du dålig aptit?
F93	Hur ofta har du under detta läsår haft ”nervös mage” (t.ex. magknip, magkramper, orolig mage, illamående, gaser, förstoppning eller diarré)?
F95	Hur ofta har du under detta läsår haft svårt att somna?
F98	Hur ofta har det hänt under detta läsår att du sovit oroligt och vaknat under natten?

Item fit

Code

RIitemfitPCM2(df.psm, 300, 32, 8)

	OutfitMSQ	InfitMSQ	OutfitZSTD	InfitZSTD
F88	0.886	0.891	-1.438	-1.2
F91	0.937	0.939	-0.624	-1.136
F93	0.904	0.901	-0.925	-1.528
F95	0.832	0.84	-2.101	-2.423
F98	0.798	0.816	-2.623	-2.55

PCA

Code

RIpcmPCA(df.psm)

PCA of Rasch model residuals

Eigenvalues
1.47
1.30
1.14
1.09
0.00

Loadings 1st contrast

Code

RIloadLoc(df.psm)

Residualkorrelationer

Code

RIresidcorr(df.psm, cutoff = 0.2)

	F88	F91	F93	F95	F98
F88
F91	-0.18
F93	-0.19	-0.08
F95	-0.27	-0.14	-0.28
F98	-0.23	-0.18	-0.26	-0.05
Note:
Relative cut-off value (highlighted in red) is 0.013, which is 0.2 above the average correlation.

Targeting

Code

RItargeting(df.psm)

Itemhierarki

Code

RIitemHierarchy(df.psm)

4.7 Invarians-/DIF

4.7.1 Kön

itemnr	item
F88	Hur ofta har du haft huvudvärk detta läsår?
F91	Hur ofta har du dålig aptit?
F93	Hur ofta har du under detta läsår haft ”nervös mage” (t.ex. magknip, magkramper, orolig mage, illamående, gaser, förstoppning eller diarré)?
F95	Hur ofta har du under detta läsår haft svårt att somna?
F98	Hur ofta har det hänt under detta läsår att du sovit oroligt och vaknat under natten?

Tabell

Code

RIdifTable(df.psm, dif.gender)

Item	2	3	Mean location	StDev	MaxDiff
F88	-0.107	0.033	-0.037	0.099	0.140
F91	0.165	0.131	0.148	0.024	0.034
F93	-0.216	-0.140	-0.178	0.054	0.076
F95	0.103	-0.130	-0.014	0.164	0.232
F98	0.056	0.105	0.080	0.035	0.050

Figur

Code

RIdifFigure(df.psm, dif.gender)

Inga problem med DIF för något item gällande kön.

4.7.2 Årskurs

itemnr	item
F88	Hur ofta har du haft huvudvärk detta läsår?
F91	Hur ofta har du dålig aptit?
F93	Hur ofta har du under detta läsår haft ”nervös mage” (t.ex. magknip, magkramper, orolig mage, illamående, gaser, förstoppning eller diarré)?
F95	Hur ofta har du under detta läsår haft svårt att somna?
F98	Hur ofta har det hänt under detta läsår att du sovit oroligt och vaknat under natten?

Tabell

Code

RIdifTable(df.psm, dif.arskurs)

Item	2	3	Mean location	StDev	MaxDiff
F88	-0.088	0.005	-0.041	0.066	0.093
F91	0.079	0.197	0.138	0.084	0.118
F93	-0.140	-0.217	-0.179	0.055	0.077
F95	-0.017	-0.007	-0.012	0.007	0.010
F98	0.166	0.022	0.094	0.102	0.145

Figur

Code

RIdifFigure(df.psm, dif.arskurs)

4.7.3 Årtal

itemnr	item
F88	Hur ofta har du haft huvudvärk detta läsår?
F91	Hur ofta har du dålig aptit?
F93	Hur ofta har du under detta läsår haft ”nervös mage” (t.ex. magknip, magkramper, orolig mage, illamående, gaser, förstoppning eller diarré)?
F95	Hur ofta har du under detta läsår haft svårt att somna?
F98	Hur ofta har det hänt under detta läsår att du sovit oroligt och vaknat under natten?

Tabell

Code

final.psm.items <- names(df.psm)
# write.csv(final.psm.items, file = "2022-09-18_PSMfinalitems.csv")
df.dif.years.psm <- df.dif.years %>%
  select(any_of(final.psm.items))

df.dif.years.psm$F91 <- recode(df.dif.years.psm$F91, "2=1;3=2;4=2", as.factor = FALSE)
df.dif.years.psm$F93 <- recode(df.dif.years.psm$F93, "4=3", as.factor = FALSE)
df.dif.years.psm$F98 <- recode(df.dif.years.psm$F98, "4=3", as.factor = FALSE)
df.dif.years.psm$F95 <- recode(df.dif.years.psm$F95, "4=3;3=2;2=1;1=0", as.factor = FALSE)

RIdifTable(df.dif.years.psm, dif.year)

Item	3	5	6	8	10	11	Mean location	StDev	MaxDiff
F88	-0.096	-0.077	-0.070	-0.043	-0.035	-0.045	-0.061	0.024	0.061
F91	0.104	0.149	0.102	0.132	0.078	0.060	0.104	0.033	0.088
F93	-0.144	-0.142	-0.125	-0.185	-0.211	-0.178	-0.164	0.032	0.085
F95	0.049	0.001	0.004	-0.009	0.038	0.002	0.014	0.024	0.058
F98	0.087	0.069	0.090	0.105	0.130	0.161	0.107	0.033	0.092

Figur 1

Code

RIdifFigure(df.dif.years.psm, dif.year)

Figur 2

Code

RIdifFigTime(df.dif.years.psm, dif.year)

Items är stabila och jämförbara över tid.

4.8 Reliabilitet - psykosomatiska besvär

Code

RItif(df.psm)

Vi har acceptabel reliabilitet för ca 68% av respondenterna.

4.9 Item-parametrar

Code

RIitemparams(df.dif.years.psm, "PsyksomItemParams.csv")

	Threshold 1	Threshold 2	Threshold 3	Threshold 4	Item location
F88	-1.09	-0.28	0.85	1.23	0.18
F91	0.01	0.70	NA	NA	0.35
F93	-0.37	0.09	0.51	NA	0.08
F95	-0.05	0.22	0.57	NA	0.25
F98	-0.08	0.25	0.86	NA	0.34

4.10 Transformeringstabell

Med mätosäkerheter angivna som Logit std.error, där värden under 0.54 motsvarar reliabilitet över PSI = 0.7.

Code

RIscoreSE(df.dif.years.psm)

Ordinal sum score	Logit score	Logit std.error
0	-2.80	NA
1	-2.02	1.00
2	-1.32	0.72
3	-0.90	0.60
4	-0.58	0.53
5	-0.32	0.49
6	-0.09	0.47
7	0.13	0.46
8	0.34	0.46
9	0.56	0.47
10	0.79	0.49
11	1.04	0.53
12	1.35	0.59
13	1.77	0.71
14	2.45	0.99
15	3.20	NA

4.11 Rasch-analys 1 - psykiska besvär

Code

items.psm <- c("F88","F91","F93","F95","F98")

df.psy <- df.omit.na %>% 
  select(!any_of(items.psm))

itemnr	item
F89	Känner du dig ledsen och deppig utan att veta varför?
F90	Händer det att du känner dig rädd utan att veta varför?
F92	Hur mycket skulle du vilja ändra på dig själv?
F94	Hur ofta tycker du att du inget duger till?
F96	Är du nöjd med ditt utseende?
F97	Känner du dig slö och olustig?
F99	Hur ofta tycker du att det är riktigt härligt att leva?

Item fit

Code

RIitemfitPCM2(df.psy, 300, 32, 8)

	OutfitMSQ	InfitMSQ	OutfitZSTD	InfitZSTD
F89	0.773	0.79	-2.851	-2.865
F90	0.958	1.007	-0.342	-0.072
F92	0.778	0.784	-2.868	-3.095
F94	0.765	0.763	-2.414	-3.208
F96	0.879	0.882	-1.293	-1.324
F97	1.031	1.021	0.336	0.253
F99	1.059	1.016	0.664	0.155

PCA

Code

RIpcmPCA(df.psy)

PCA of Rasch model residuals

Eigenvalues
1.60
1.28
1.18
1.09
0.98

Residualkorrelationer

Code

RIresidcorr(df.psy, cutoff = 0.2)

	F89	F90	F92	F94	F96	F97	F99
F89
F90	0.06
F92	-0.17	-0.18
F94	-0.11	-0.15	-0.08
F96	-0.34	-0.26	0.12	-0.05
F97	-0.14	-0.12	-0.24	-0.16	-0.2
F99	-0.12	-0.18	-0.18	-0.22	-0.11	-0.16
Note:
Relative cut-off value (highlighted in red) is 0.058, which is 0.2 above the average correlation.

Targeting

Code

RItargeting(df.psy)

Itemhierarki

Code

RIitemHierarchy(df.psy)

Svarskategorier

Code

plot(mirt.rasch, type="trace")

Barplots

Code

for (i in 1:ncol(df.psy)) {
    barplot(table(df.psy[, i]), col = "#8dc8c7", main = names(df.psy[i]), 
      ylab = "Number of responses", xlab = (itemlabels %>% 
                                              filter(itemnr %in% names(df.psy)) 
                                            %>% .[i,2]))
  }

Items

Code

RIlistitems(df.psy)

itemnr	item
F89	Känner du dig ledsen och deppig utan att veta varför?
F90	Händer det att du känner dig rädd utan att veta varför?
F92	Hur mycket skulle du vilja ändra på dig själv?
F94	Hur ofta tycker du att du inget duger till?
F96	Är du nöjd med ditt utseende?
F97	Känner du dig slö och olustig?
F99	Hur ofta tycker du att det är riktigt härligt att leva?

Många items har problem med oordnade svarskategorier.

Items 92 och 96 har en för stor residualkorrelation, men vi avvaktar med åtgärd tills svarskategorierna är åtgärdade.

4.12 Omkodning av svarskategorier

Vi behöver åtgärda flera items:

• F89, 90, 96: 1+2 och 3+4
• F92, 99: 3+4
• F97: 0+1

(flera varianter testades, och ovanstående fungerade minst dåligt)

Code

df.psy$F89 <- recode(df.psy$F89, "2=1;3=2;4=2", as.factor = FALSE)
df.psy$F90 <- recode(df.psy$F90, "2=1;3=2;4=2", as.factor = FALSE)
df.psy$F96 <- recode(df.psy$F96, "2=1;3=2;4=2", as.factor = FALSE)
df.psy$F97 <- recode(df.psy$F97, "1=0;2=1;3=2;4=3", as.factor = FALSE)

df.psy$F92 <- recode(df.psy$F92, "4=3", as.factor = FALSE)
df.psy$F99 <- recode(df.psy$F99, "4=3", as.factor = FALSE)

Svarskategorier

Code

plot(mirt.rasch, type="trace")

Itemhierarki (med konfidensintervall runt svarströsklar)

Code

RIitemHierarchy(df.psy)

4.13 Rasch-analys 2 - psykiska besvär

Efter att ha åtgärdat svarskategorierna.

itemnr	item
F89	Känner du dig ledsen och deppig utan att veta varför?
F90	Händer det att du känner dig rädd utan att veta varför?
F92	Hur mycket skulle du vilja ändra på dig själv?
F94	Hur ofta tycker du att du inget duger till?
F96	Är du nöjd med ditt utseende?
F97	Känner du dig slö och olustig?
F99	Hur ofta tycker du att det är riktigt härligt att leva?

Item fit

Code

RIitemfitPCM2(df.psy, 300, 32, 8)

	OutfitMSQ	InfitMSQ	OutfitZSTD	InfitZSTD
F89	0.777	0.79	-3.002	-3.035
F90	0.912	0.951	-0.68	-0.477
F92	0.832	0.841	-1.909	-2.269
F94	0.774	0.792	-2.438	-2.786
F96	0.845	0.845	-2.073	-2.126
F97	1.071	1.024	0.707	0.404
F99	1.076	1.045	0.762	0.487

PCA

Code

RIpcmPCA(df.psy)

PCA of Rasch model residuals

Eigenvalues
1.52
1.33
1.21
1.13
0.91

Residualkorrelationer

Code

RIresidcorr(df.psy, cutoff = 0.2)

	F89	F90	F92	F94	F96	F97	F99
F89
F90	0.12
F92	-0.13	-0.13
F94	-0.11	-0.13	-0.17
F96	-0.2	-0.18	0.12	-0.07
F97	-0.11	-0.11	-0.24	-0.22	-0.17
F99	-0.11	-0.16	-0.24	-0.31	-0.12	-0.16
Note:
Relative cut-off value (highlighted in red) is 0.065, which is 0.2 above the average correlation.

Targeting

Code

RItargeting(df.psy)

Itemhierarki

Code

RIitemHierarchy(df.psy)

Vi har två residualkorrelationer ca 0.25 över medelvärdet för samtliga residualkorrelationer.

• F96 tas bort (F92 behålls pga bättre targeting)
• F90 tas bort (F89 behålls pga bättre targeting)

4.13.1 Reliabilitet

Även om vi inte undersökt DIF, vilket bör göras innan reliabiliteten bedöms, så har de påtagliga problemen med svarskategorierna medfört att vi kan ha såpass låg reliabilitet att det inte är meningsfullt med DIF-analysen.

Code

df.psy$F96 <- NULL
df.psy$F90 <- NULL

RItif(df.psy)

Eftersom många items uppvisade problem med svarskategorierna blir reliabiliteten låg.

Vi provar att slå samman kvarvarande items från psykiska besvär med de items som utgör psykosomatiska besvär.

Code

df.psf <- cbind(df.psm,df.psy)

df.psf <- df.psf |> 
  relocate(F92, .after = "F91") |> 
  relocate(F97, .before = "F98") |> 
  relocate(F89, .after = "F88") |>
  relocate(F94, .before = "F95")

# psf.items <- names(df.psf)
# 
# df.psf <- df.omit.na |> 
#   select(any_of(psf.items))

4.14 Rasch-analys 1 psykiska/psykosomatiska besvär

itemnr	item
F88	Hur ofta har du haft huvudvärk detta läsår?
F89	Känner du dig ledsen och deppig utan att veta varför?
F91	Hur ofta har du dålig aptit?
F92	Hur mycket skulle du vilja ändra på dig själv?
F93	Hur ofta har du under detta läsår haft ”nervös mage” (t.ex. magknip, magkramper, orolig mage, illamående, gaser, förstoppning eller diarré)?
F94	Hur ofta tycker du att du inget duger till?
F95	Hur ofta har du under detta läsår haft svårt att somna?
F97	Känner du dig slö och olustig?
F98	Hur ofta har det hänt under detta läsår att du sovit oroligt och vaknat under natten?
F99	Hur ofta tycker du att det är riktigt härligt att leva?

Item fit

Code

RIitemfitPCM2(df.psf, 300, 32, 8)

	OutfitMSQ	InfitMSQ	OutfitZSTD	InfitZSTD
F88	1.078	1.058	1.263	0.821
F89	0.719	0.728	-4.626	-4.358
F91	0.963	0.961	-0.507	-0.437
F92	0.86	0.863	-1.722	-2.189
F93	1.014	0.993	0.468	-0.336
F94	0.904	0.887	-1.003	-1.456
F95	1.006	0.985	-0.279	0.111
F97	0.967	0.954	-0.201	-0.423
F98	0.95	0.957	-0.996	-0.749
F99	1.064	1.01	0.696	-0.026

PCA

Code

RIpcmPCA(df.psf)

PCA of Rasch model residuals

Eigenvalues
1.78
1.26
1.13
1.08
1.02

Loadings 1st contrast

Code

RIloadLoc(df.psf)

Residualkorrelationer

Code

RIresidcorr(df.psf, cutoff = 0.2)

	F88	F89	F91	F92	F93	F94	F95	F97	F98	F99
F88
F89	-0.1
F91	-0.07	-0.03
F92	-0.22	0.05	-0.06
F93	-0.03	-0.07	0	-0.1
F94	-0.23	0.08	-0.16	0.11	-0.18
F95	-0.08	-0.14	-0.04	-0.21	-0.11	-0.26
F97	-0.14	0.01	-0.13	-0.05	-0.18	0.02	-0.13
F98	-0.05	-0.09	-0.08	-0.2	-0.1	-0.21	0.11	-0.17
F99	-0.23	0.04	-0.12	0.01	-0.2	-0.02	-0.14	0.01	-0.19
Note:
Relative cut-off value (highlighted in red) is 0.109, which is 0.2 above the average correlation.

Targeting

Code

RItargeting(df.psf)

Itemhierarki

Code

RIitemHierarchy(df.psf)

Svarskategorier

Code

plot(mirt.rasch, type="trace")

ICC plot

Code

df.psf$F89 <- NULL
df.psf$F94 <- recode(df.psf$F94, "4=3;3=2;2=1;1=0", as.factor = FALSE)
df.erm <- PCM(df.psf)
plotICC(df.erm, item.subset = "F94")

Konfidensintervall för svarströsklar

Code

RIitemHierarchy(df.psf)

F94 har fortfarande problem med svarskategorierna. Vi slår samman de två mittersta.

ICC plot

Code

df.psf$F94 <- recode(df.psf$F94, "3=2;2=1", as.factor = FALSE)
df.erm <- PCM(df.psf)
plotICC(df.erm, item.subset = "F94")

Konfidensintervall för svarströsklar

Code

RIitemHierarchy(df.psf)

4.15 Rasch-analys 2 psykiska/psykosomatiska besvär

itemnr	item
F88	Hur ofta har du haft huvudvärk detta läsår?
F91	Hur ofta har du dålig aptit?
F92	Hur mycket skulle du vilja ändra på dig själv?
F93	Hur ofta har du under detta läsår haft ”nervös mage” (t.ex. magknip, magkramper, orolig mage, illamående, gaser, förstoppning eller diarré)?
F94	Hur ofta tycker du att du inget duger till?
F95	Hur ofta har du under detta läsår haft svårt att somna?
F97	Känner du dig slö och olustig?
F98	Hur ofta har det hänt under detta läsår att du sovit oroligt och vaknat under natten?
F99	Hur ofta tycker du att det är riktigt härligt att leva?

Item fit

Code

RIitemfitPCM2(df.psf, 300, 32, 8)

	OutfitMSQ	InfitMSQ	OutfitZSTD	InfitZSTD
F88	1.009	0.997	0.038	-0.06
F91	0.929	0.936	-0.919	-1.045
F92	0.871	0.874	-1.736	-1.727
F93	0.961	0.95	-0.4	-0.571
F94	0.835	0.839	-1.932	-2.329
F95	0.922	0.918	-0.992	-1.229
F97	0.959	0.948	-0.314	-0.524
F98	0.895	0.908	-0.992	-1.489
F99	1.058	1.009	0.972	0.213

PCA

Code

RIpcmPCA(df.psf)

PCA of Rasch model residuals

Eigenvalues
1.69
1.29
1.16
1.07
1.01

Loadings 1st contrast

Code

RIloadLoc(df.psf)

Residualkorrelationer

Code

RIresidcorr(df.psf, cutoff = 0.2)

	F88	F91	F92	F93	F94	F95	F97	F98	F99
F88
F91	-0.09
F92	-0.22	-0.05
F93	-0.06	-0.02	-0.1
F94	-0.14	-0.09	0.16	-0.11
F95	-0.12	-0.07	-0.21	-0.15	-0.16
F97	-0.15	-0.13	-0.03	-0.18	0.09	-0.15
F98	-0.08	-0.1	-0.2	-0.13	-0.12	0.07	-0.18
F99	-0.24	-0.12	0.03	-0.2	0.05	-0.16	0.01	-0.2
Note:
Relative cut-off value (highlighted in red) is 0.102, which is 0.2 above the average correlation.

Targeting

Code

RItargeting(df.psf)

Itemhierarki

Code

RIitemHierarchy(df.psf)

Svarskategorier

Code

plot(mirt.rasch, type="trace")

F94 korrelerar för starkt med F92. Vi tar bort F94 p.g.a. sämre targeting och ZSTD.

Code

df.psf$F94 <- NULL

4.16 Rasch-analys 3 psykiska/psykosomatiska besvär

itemnr	item
F88	Hur ofta har du haft huvudvärk detta läsår?
F91	Hur ofta har du dålig aptit?
F92	Hur mycket skulle du vilja ändra på dig själv?
F93	Hur ofta har du under detta läsår haft ”nervös mage” (t.ex. magknip, magkramper, orolig mage, illamående, gaser, förstoppning eller diarré)?
F95	Hur ofta har du under detta läsår haft svårt att somna?
F97	Känner du dig slö och olustig?
F98	Hur ofta har det hänt under detta läsår att du sovit oroligt och vaknat under natten?
F99	Hur ofta tycker du att det är riktigt härligt att leva?

Item fit

Code

RIitemfitPCM2(df.psf, 300, 32, 8)

	OutfitMSQ	InfitMSQ	OutfitZSTD	InfitZSTD
F88	0.962	0.958	-0.554	-0.566
F91	0.902	0.915	-1.238	-1.462
F92	0.886	0.885	-1.667	-1.734
F93	0.923	0.922	-0.612	-1.136
F95	0.882	0.881	-1.355	-1.819
F97	0.972	0.947	-0.68	-0.84
F98	0.863	0.878	-1.736	-1.777
F99	1.055	1.004	0.646	0.296

PCA

Code

RIpcmPCA(df.psf)

PCA of Rasch model residuals

Eigenvalues
1.55
1.31
1.17
1.05
1.00

Loadings 1st contrast

Code

RIloadLoc(df.psf)

Residualkorrelationer

Code

RIresidcorr(df.psf, cutoff = 0.2)

	F88	F91	F92	F93	F95	F97	F98	F99
F88
F91	-0.1
F92	-0.21	-0.04
F93	-0.07	-0.02	-0.08
F95	-0.14	-0.08	-0.2	-0.16
F97	-0.14	-0.13	0	-0.17	-0.14
F98	-0.09	-0.11	-0.18	-0.14	0.06	-0.17
F99	-0.24	-0.12	0.04	-0.2	-0.16	0.03	-0.2
Note:
Relative cut-off value (highlighted in red) is 0.087, which is 0.2 above the average correlation.

Targeting

Code

RItargeting(df.psf)

Itemhierarki

Code

RIitemHierarchy(df.psf)

Detta ser ut att fungera bra. Vi övergår till DIF-analys

4.17 DIF/invarians

4.17.1 Kön

itemnr	item
F88	Hur ofta har du haft huvudvärk detta läsår?
F91	Hur ofta har du dålig aptit?
F92	Hur mycket skulle du vilja ändra på dig själv?
F93	Hur ofta har du under detta läsår haft ”nervös mage” (t.ex. magknip, magkramper, orolig mage, illamående, gaser, förstoppning eller diarré)?
F95	Hur ofta har du under detta läsår haft svårt att somna?
F97	Känner du dig slö och olustig?
F98	Hur ofta har det hänt under detta läsår att du sovit oroligt och vaknat under natten?
F99	Hur ofta tycker du att det är riktigt härligt att leva?

Tabell

Code

#dif.gender<-recode(dif.gender,"'Pojke'=1;'Flicka'=2",as.factor=FALSE)
RIdifTable(df.psf, dif.gender)

Item	2	3	Mean location	StDev	MaxDiff
F88	-0.204	-0.012	-0.108	0.136	0.193
F91	0.056	0.086	0.071	0.021	0.030
F92	-0.623	-0.443	-0.533	0.127	0.180
F93	-0.308	-0.180	-0.244	0.091	0.128
F95	-0.002	-0.171	-0.086	0.119	0.169
F97	0.779	0.628	0.704	0.106	0.151
F98	-0.048	0.056	0.004	0.074	0.104
F99	0.350	0.035	0.193	0.223	0.315

Figur

Code

RIdifFigure(df.psf, dif.gender)

Svarskategorier

Code

RIdifFigThresh(df.psf, dif.gender)

F99 är enda item som går strax över 0.3 logits. Behöver ej åtgärdas.

4.17.2 Årskurs

itemnr	item
F88	Hur ofta har du haft huvudvärk detta läsår?
F91	Hur ofta har du dålig aptit?
F92	Hur mycket skulle du vilja ändra på dig själv?
F93	Hur ofta har du under detta läsår haft ”nervös mage” (t.ex. magknip, magkramper, orolig mage, illamående, gaser, förstoppning eller diarré)?
F95	Hur ofta har du under detta läsår haft svårt att somna?
F97	Känner du dig slö och olustig?
F98	Hur ofta har det hänt under detta läsår att du sovit oroligt och vaknat under natten?
F99	Hur ofta tycker du att det är riktigt härligt att leva?

Tabell

Code

RIdifTable(df.psf, dif.arskurs)

Item	2	3	Mean location	StDev	MaxDiff
F88	-0.157	-0.077	-0.117	0.056	0.079
F91	0.004	0.111	0.057	0.076	0.107
F92	-0.580	-0.517	-0.548	0.044	0.062
F93	-0.208	-0.286	-0.247	0.055	0.078
F95	-0.090	-0.085	-0.088	0.003	0.005
F97	0.744	0.699	0.722	0.032	0.045
F98	0.086	-0.058	0.014	0.102	0.144
F99	0.200	0.214	0.207	0.010	0.014

Figur

Code

RIdifFigure(df.psf, dif.arskurs)

Svarskategorier

Code

RIdifFigThresh(df.psf, dif.arskurs)

Inga skillnader mellan årskurser.

4.17.3 Årtal

itemnr	item
F88	Hur ofta har du haft huvudvärk detta läsår?
F91	Hur ofta har du dålig aptit?
F92	Hur mycket skulle du vilja ändra på dig själv?
F93	Hur ofta har du under detta läsår haft ”nervös mage” (t.ex. magknip, magkramper, orolig mage, illamående, gaser, förstoppning eller diarré)?
F95	Hur ofta har du under detta läsår haft svårt att somna?
F97	Känner du dig slö och olustig?
F98	Hur ofta har det hänt under detta läsår att du sovit oroligt och vaknat under natten?
F99	Hur ofta tycker du att det är riktigt härligt att leva?

Tabell

Code

final.items.psf <- names(df.psf)
#write.csv(final.items, file = "2022-09-16 IFoptimalItems.csv")
df.dif.years <- df.dif.years %>% 
  select(any_of(final.items.psf))

df.dif.years$F91 <- recode(df.dif.years$F91, "2=1;3=2;4=2", as.factor = FALSE)
df.dif.years$F92 <- recode(df.dif.years$F92, "4=3", as.factor = FALSE)
df.dif.years$F93 <- recode(df.dif.years$F93, "4=3", as.factor = FALSE)
df.dif.years$F98 <- recode(df.dif.years$F98, "4=3", as.factor = FALSE)
df.dif.years$F99 <- recode(df.dif.years$F99, "4=3", as.factor = FALSE)
df.dif.years$F95 <- recode(df.dif.years$F95, "4=3;3=2;2=1;1=0", as.factor = FALSE)
df.dif.years$F97 <- recode(df.dif.years$F97, "4=3;3=2;2=1;1=0", as.factor = FALSE)

RIdifTable(df.dif.years, dif.year)

Item	3	6	7	8	11	12	14	15	Mean location	StDev	MaxDiff
F88	-0.179	-0.177	-0.180	-0.161	-0.111	-0.148	-0.112	-0.141	-0.151	0.028	0.070
F91	0.017	0.046	0.039	0.007	0.060	0.022	-0.003	-0.036	0.019	0.030	0.096
F92	-0.502	-0.504	-0.494	-0.518	-0.542	-0.508	-0.495	-0.503	-0.508	0.016	0.048
F93	-0.221	-0.241	-0.233	-0.210	-0.247	-0.277	-0.278	-0.266	-0.247	0.025	0.068
F95	-0.038	-0.121	-0.082	-0.088	-0.088	-0.102	-0.042	-0.092	-0.082	0.028	0.084
F97	0.777	0.811	0.781	0.744	0.719	0.769	0.693	0.741	0.754	0.038	0.118
F98	-0.003	-0.045	-0.029	-0.007	0.002	0.021	0.045	0.061	0.006	0.035	0.106
F99	0.149	0.232	0.198	0.234	0.208	0.223	0.193	0.234	0.209	0.029	0.086

Figur items

Code

RIdifFigure(df.dif.years, dif.year)

Över tid

Code

RIdifFigTime(df.dif.years, dif.year)

Svarskategorier

Code

RIdifFigThresh(df.dif.years, dif.year)

Items fungerar stabilt över tid.

4.17.4 Interaktion årskurs och kön

itemnr	item
F88	Hur ofta har du haft huvudvärk detta läsår?
F91	Hur ofta har du dålig aptit?
F92	Hur mycket skulle du vilja ändra på dig själv?
F93	Hur ofta har du under detta läsår haft ”nervös mage” (t.ex. magknip, magkramper, orolig mage, illamående, gaser, förstoppning eller diarré)?
F95	Hur ofta har du under detta läsår haft svårt att somna?
F97	Känner du dig slö och olustig?
F98	Hur ofta har det hänt under detta läsår att du sovit oroligt och vaknat under natten?
F99	Hur ofta tycker du att det är riktigt härligt att leva?

Code

dfin <- df.psf

df.tree <- data.frame(matrix(ncol = 0, nrow = nrow(dfin)))
df.tree$difdata <- as.matrix(dfin)
df.tree$dif.gender <- dif.gender
df.tree$dif.arskurs <- dif.arskurs

pctree.out <- pctree(difdata ~ dif.gender + dif.arskurs, data = df.tree)
plot(pctree.out)

Code

cutoff <- 0.5
diffig <- itempar(pctree.out) %>%
  as.data.frame() %>%
  t() %>%
  as.data.frame() %>%
  mutate(
    `Mean location` = rowMeans(.),
    StDev = rowSds(as.matrix(.))
  ) %>%
  rowwise() %>%
  mutate(MaxDiff = (max(c_across(c(1:(ncol(.) - 2))))) -
    min(c_across(c(1:(ncol(.) - 2))))) %>%
  ungroup() %>%
  mutate(across(where(is.numeric), round, 3)) %>%
  rownames_to_column(var = "Item") %>%
  mutate(Item = names(dfin)) %>%
  relocate(MaxDiff, .after = last_col())

  formattable(diffig,
    list(MaxDiff = formatter("span",
      style = ~ style(color = ifelse(MaxDiff < -cutoff,
        "red", ifelse(MaxDiff > cutoff, "red", "black")
      ))
    ),
    formattable::area(col = 2:6) ~ color_tile(RISEprimYellowLight, RISEprimGreen)
    ),
    table.attr = "class=\"table table-striped\" style=\"font-size: 15px; font-family: Lato\""
  )

Item	3	4	6	7	Mean location	StDev	MaxDiff
F88	-0.235	-0.184	-0.071	0.036	-0.113	0.121	0.271
F91	-0.029	0.114	0.046	0.119	0.062	0.070	0.149
F92	-0.692	-0.576	-0.439	-0.446	-0.538	0.120	0.253
F93	-0.241	-0.355	-0.167	-0.190	-0.238	0.083	0.188
F95	0.019	-0.018	-0.183	-0.160	-0.085	0.101	0.202
F97	0.820	0.753	0.627	0.628	0.707	0.096	0.193
F98	0.075	-0.133	0.092	0.028	0.016	0.103	0.225
F99	0.283	0.399	0.095	-0.015	0.190	0.186	0.414

4.17.5 Item fit över tid

DIF-analyserna visar på stabilitet i item locations över tid. Det kan också vara av intresse att se hur item fit ser ut över tid (2006-2020).

Code

final.items <- read_csv("PSFitemnr.csv") %>% 
  pull(itemnr)

df.fit <- df.all %>%
  select(final.items,ar) %>% 
  na.omit()

df.fit$F91 <- recode(df.fit$F91, "2=1;3=2;4=2", as.factor = FALSE)
df.fit$F92 <- recode(df.fit$F92, "4=3", as.factor = FALSE)
df.fit$F93 <- recode(df.fit$F93, "4=3", as.factor = FALSE)
df.fit$F98 <- recode(df.fit$F98, "4=3", as.factor = FALSE)
df.fit$F99 <- recode(df.fit$F99, "4=3", as.factor = FALSE)
df.fit$F95 <- recode(df.fit$F95, "4=3;3=2;2=1;1=0", as.factor = FALSE)
df.fit$F97 <- recode(df.fit$F97, "4=3;3=2;2=1;1=0", as.factor = FALSE)

Code

library(furrr)
plan(multisession, workers = 8)

itemfit.time <- df.fit %>% 
  split(.$ar) %>% 
  future_map(~ RIitemfitPCM(subset(., select = -ar), table = FALSE), 
             .options = furrr_options(seed = TRUE))

årtal <- c("2006", "2008", "2010", "2012", "2014", "2016", "2018", "2020")
itemfit.all <- data.frame(matrix(ncol = 4, nrow = 0))
colnames(itemfit.all) <- c("Item", "År", 'InfitMSQ','OutfitMSQ')

for (i in årtal) {
  itemfit.add <- itemfit.time[[i]] %>% 
    rownames_to_column("Item") %>% 
    add_column(År = i)
  
  itemfit.all <- rbind(itemfit.all, itemfit.add)
}

Code

itemfit.all %>%
  group_by(Item) %>%
  summarise(
    mean_Out = mean(OutfitMSQ),
    sd_Out = sd(OutfitMSQ),
    max_Out = max(OutfitMSQ),
    min_Out = min(OutfitMSQ),
    mean_In = mean(InfitMSQ),
    sd_In = sd(InfitMSQ),
    max_In = max(InfitMSQ),
    min_In = min(InfitMSQ)
  ) %>%
  mutate(
    diffmax_OutfitMSQ = max_Out - min_Out,
    diffmax_InfitMSQ = max_In - min_In,
    item = Item
  ) %>%
  mutate(across(where(is.numeric), ~ round(.x, 2))) %>%
  mutate(across(starts_with(c("max", "min")), ~ cell_spec(.x, color = ifelse(.x < msq_min, "red",
    ifelse(.x > msq_max, "red", "black")
  )))) %>%
  kbl(booktabs = T, escape = F,
      col.names = c("Item","Mean","SD","Max","Min","Mean","SD","Max","Min","OutfitMSQ","InfitMSQ","Item")) %>%
  # bootstrap options are for HTML output
  kable_styling(
    bootstrap_options = c("striped", "hover"),
    position = "left",
    full_width = F,
    font_size = 15,
    fixed_thead = T
  ) %>%
  column_spec(c(1, 12), bold = T) %>%
  kable_classic(html_font = "Lato") %>%
  # latex_options are for PDF output
  kable_styling(latex_options = c("striped", "scale_down")) %>%
  add_header_above(c(" " = 1, 
                     "OutfitMSQ" = 4, 
                     "InfitMSQ" = 4, 
                     "Max differences" = 2,
                     " " = 1))

	OutfitMSQ				InfitMSQ				Max differences
Item	Mean	SD	Max	Min	Mean	SD	Max	Min	OutfitMSQ	InfitMSQ	Item
F88	0.97	0.02	1	0.96	0.96	0.01	0.99	0.95	0.04	0.04	F88
F91	0.88	0.01	0.91	0.86	0.89	0.01	0.92	0.87	0.05	0.04	F91
F92	0.91	0.02	0.95	0.89	0.91	0.02	0.94	0.89	0.06	0.05	F92
F93	0.93	0.01	0.95	0.91	0.93	0.01	0.94	0.91	0.04	0.03	F93
F95	0.88	0.03	0.94	0.85	0.89	0.02	0.93	0.86	0.09	0.07	F95
F97	0.95	0.03	0.99	0.9	0.94	0.02	0.96	0.9	0.09	0.06	F97
F98	0.86	0.01	0.86	0.84	0.87	0.01	0.88	0.86	0.02	0.02	F98
F99	1.04	0.02	1.06	1.01	1.00	0.01	1.01	0.98	0.05	0.03	F99

Inget item går utanför gränsvärden 0.7-1.3.

4.18 Reliabilitet - psykiska/psykosomatiska besvär

Code

RItif(df.psf)

Bättre reliabilitet än vad någon av delskalorna uppnår separat.

4.19 Person fit

Code

RIpfit(df.psf)

4.20 Item-parametrar

itemnr	item
F88	Hur ofta har du haft huvudvärk detta läsår?
F91	Hur ofta har du dålig aptit?
F92	Hur mycket skulle du vilja ändra på dig själv?
F93	Hur ofta har du under detta läsår haft ”nervös mage” (t.ex. magknip, magkramper, orolig mage, illamående, gaser, förstoppning eller diarré)?
F95	Hur ofta har du under detta läsår haft svårt att somna?
F97	Känner du dig slö och olustig?
F98	Hur ofta har det hänt under detta läsår att du sovit oroligt och vaknat under natten?
F99	Hur ofta tycker du att det är riktigt härligt att leva?

Code

RIitemparams(df.dif.years)

	Threshold 1	Threshold 2	Threshold 3	Threshold 4	Item location
F88	-1.09	-0.37	0.74	1.01	0.07
F91	-0.06	0.55	NA	NA	0.25
F92	-1.36	0.06	0.46	NA	-0.28
F93	-0.41	-0.01	0.36	NA	-0.02
F95	-0.10	0.13	0.41	NA	0.14
F97	0.36	1.06	1.52	NA	0.98
F98	-0.14	0.14	0.68	NA	0.23
F99	-0.32	0.60	1.03	NA	0.44

Code

itemlabels %>% 
  filter(itemnr %in% names(df.dif.years)) %>% 
  write_csv("PSFitemnr.csv")

Code

RIscoreSE(df.dif.years)

Ordinal sum score	Logit score	Logit std.error
0	-3.38	NA
1	-2.59	1.01
2	-1.88	0.72
3	-1.45	0.60
4	-1.14	0.52
5	-0.89	0.47
6	-0.68	0.44
7	-0.50	0.41
8	-0.34	0.40
9	-0.19	0.38
10	-0.04	0.38
11	0.10	0.37
12	0.24	0.37
13	0.37	0.37
14	0.51	0.38
15	0.66	0.38
16	0.81	0.39
17	0.97	0.41
18	1.14	0.43
19	1.34	0.46
20	1.58	0.51
21	1.87	0.58
22	2.27	0.70
23	2.94	0.99
24	3.69	NA

4.21 Exempel med single-item vs index

Data från enbart 2020.

OBS! Viktigt att komma ihåg att höga värden = hög risk, vilket även gäller svar på enskilda frågor. En positiv fråga som F99 har alltså omvända svarskategorier (3 och 4 sammanslagna):

• 3-4 - Sällan / Någon enstaka gång
• 2 - Ibland
• 1 - Ganska ofta
• 0 - Väldigt ofta

Code

#---- create df with only 2020 respondents as example----
df.if <- df %>% 
  filter(ar == 2020) %>% 
  select(any_of(names(df.psf)),Kön) %>% 
  na.omit()

F99o <- df.if$F99

# rcat2 <- c("F95","F97")
# rcat3 <- c("F92","F93","F99")
# 
# for (i in rcat2) {
#   df.if[[i]]<-recode(df.if[[i]],"1=0;2=1;3=2;4=3",as.factor=FALSE)
# }
# for (i in rcat3) {
#   df.if[[i]]<-recode(df.if[[i]],"4=3",as.factor=FALSE)
# }
df.if$F91 <- recode(df.if$F91, "2=1;3=2;4=2", as.factor = FALSE)
df.if$F92 <- recode(df.if$F92, "4=3", as.factor = FALSE)
df.if$F93 <- recode(df.if$F93, "4=3", as.factor = FALSE)
df.if$F95 <- recode(df.if$F95, "4=3;3=2;2=1;1=0", as.factor = FALSE)
df.if$F97 <- recode(df.if$F97, "4=3;3=2;2=1;1=0", as.factor = FALSE)
df.if$F98 <- recode(df.if$F98, "4=3", as.factor = FALSE)
df.if$F99 <- recode(df.if$F99, "4=3", as.factor = FALSE)
df.if.gender <- df.if$Kön
df.if$Kön <- NULL
#---- estimate person locations----
library(catR)

# create matrix with item parameters
df.erm <- PCM(df.psf)
item.estimates <- eRm::thresholds(df.erm)
item_difficulty <- as.data.frame(item.estimates[["threshtable"]][["1"]])
item_difficulty$Location <- NULL
items <- item_difficulty %>% 
  mutate_if(is.character, as.numeric) %>% 
  as.matrix()

# loop through all participants to calculate theta for each one
thetaEstScores <- c()
#p1 <- c()
for (i in 1:nrow(df.if)){
  p1 <- as.numeric(as.vector(df.if[i,]))
  ptheta <- thetaEst(items, p1, model = "PCM", method = "ML") # several method options available for theta calculation
  thetaEstScores <- c(thetaEstScores,ptheta)
}

# # loop through all participants theta scores to calculate SE for each one
# thetaEstSE <- c()
# for (i in 1:length(thetaEstScores)) {
#   p1theta <- thetaEstScores[i]
#   p1 <- as.numeric(as.vector(df.if[i,]))
#   pthetaSE <- semTheta(thEst = p1theta, it = items, x = p1, model = "PCM", method = "ML") # several method options available for theta calculation
#   thetaEstSE <- c(thetaEstSE,pthetaSE)
# }

df.if$SumScores <- thetaEstScores
df.if$Kön <- df.if.gender
df.if$F99o <- F99o
#write_parquet(df.if, sink = "PSF2020.parquet")

#---- barplot for single item----
barplot(table(df.if$F99), col = "#8dc8c7", main = "Svarsfördelning för item F99")

Code

#---- visualize single item vs index score----
# some code from https://dallasnova.rbind.io/post/efficient-data-visualization-with-faded-raincloud-plots-delete-boxplot/
library(ggdist) # for shadeable density slabs
library(gghalves) # for half-half geoms
library(ggpp) # for position_dodge2nudge
library(colorspace) # for lightening color palettes
library(extrafont) # for Lato, RISE font
library(stringr)
# Setting colorblind-friendly palette
cbPalette <-c("#999999","#E69F00", "#56B4E9","#009E73",
              "#F0E442", "#0072B2", "#D55E00","#CC79A7")

cbPalette <- lighten(cbPalette, amount = 0.6, space = "HLS")

#---- plot for index value distributions divided by gender----

ggplot(df.if, aes(x = Kön, y = SumScores)) +
  stat_slab(side = "right", show.legend = F,
            scale = 0.6, # defines the height that a slab can reach
            position = position_dodge(width=.6), # distance between elements for dodging
            aes(fill_ramp = stat(level), fill=Kön), 
            .width = c(.50,.95,1)) +  # set shading
  stat_summary(fun.data = "mean_cl_normal",show.legend = F, size = .4,
               position = position_dodge2nudge(x=.05,width = .8)) +
# styling
  scale_fill_ramp_discrete(from='black', aesthetics = "fill_ramp")+ # set ramping color
  guides( # change name and display of legend elements
          color="none") + # suppresses color legend item) 
  scale_colour_manual(values = cbPalette, aesthetics = c("colour","fill"))+
  theme(plot.caption = element_text(face = "italic")) +
  #xlab(str_wrap("Svarskategorier för item F99", width = 40)) +
  ylab("Indexvärde") +
  ggtitle("Svarsfördelning för psykisk och psykosomatisk hälsa") +
  #labs(caption = "Punkterna visar medelvärden för indexvärde, fälten motsvarar 50% och 95% av distributionen.",
  #     subtitle = "Item F99 - Hur ofta tycker du att det är riktigt härligt att leva?") +
  coord_flip()

Code

#---- plots for index value vs single item response ----
ggplot(df.if, aes(x = as.factor(F99), y = SumScores)) +
  #geom_violin()
  # density distribution slab
  #geom_boxplot(size = 0.2) +
  stat_slab(side = "right", show.legend = F,
            scale = 0.6, # defines the height that a slab can reach
            position = position_dodge(width=.6), # distance between elements for dodging
            aes(fill_ramp = stat(level), fill=as.factor(F99)), 
            .width = c(.50,.95,1)) +  # set shading
  stat_summary(fun.data = "mean_cl_normal", show.legend = F, size = .4,
               position = position_dodge2nudge(x=.05,width = .8)) +
# styling
  scale_fill_ramp_discrete(from='black', aesthetics = "fill_ramp")+ # set ramping color
  guides( # change name and display of legend elements
          color="none") + # suppresses color legend item) 
  scale_colour_manual(values = cbPalette, aesthetics = c("colour","fill"))+
  theme(plot.caption = element_text(face = "italic")) +
  xlab(str_wrap("Svarskategorier för item F99", width = 40)) +
  ylab("Indexvärde") +
  ggtitle("Jämförelse av ett item och indexvärde") +
  labs(caption = "Punkterna visar medelvärden för indexvärde, fälten motsvarar 50% och 95% av distributionen.",
       subtitle = "Item F99 - Hur ofta tycker du att det är riktigt härligt att leva?") +
  coord_flip()

Code

ggplot(df.if, aes(x = as.factor(F99), y = SumScores)) +
  #geom_violin()
  # density distribution slab
  #geom_boxplot(size = 0.2) +
  stat_slab(side = "right", show.legend = T,
            scale = 0.6, # defines the height that a slab can reach
            position = position_dodge(width=.6), # distance between elements for dodging
            aes(fill_ramp = stat(level), fill=Kön), 
            .width = c(.50,.95,1)) +  # set shading
  stat_summary(fun.data = "mean_cl_normal",show.legend = F, size = .4,
               position = position_dodge2nudge(x=.05,width = .8)) +
# styling
  scale_fill_ramp_discrete(from='black', aesthetics = "fill_ramp")+ # set ramping color
  guides( # change name and display of legend elements
          color="none") + # suppresses color legend item) 
  scale_colour_manual(values = cbPalette, aesthetics = c("colour","fill"))+
  theme(plot.caption = element_text(face = "italic")) +
  xlab(str_wrap("Svarskategorier för item F99", width = 40)) +
  ylab("Indexvärde") +
  ggtitle("Jämförelse av ett item och indexvärde") +
  labs(caption = "Punkterna visar medelvärden för indexvärde, fälten motsvarar 50% och 95% av distributionen.",
       subtitle = "Item F99 - Hur ofta tycker du att det är riktigt härligt att leva?") +
  coord_flip()

Code

# stat_dotsinterval(scale = 0.5, quantiles = 100, position = "dodge") +

4.21.1 Distribution med ej omkodade svarskategorier

Code

barplot(table(df.if$F99o), col = "#8dc8c7", main = "Svarsfördelning för item F99 innan recode")

Code

#plotICC(df.ermo, item.subset = "F99")

ggplot(df.if, aes(x = as.factor(F99o), y = SumScores)) +
  #geom_violin()
  # density distribution slab
  #geom_boxplot(size = 0.2) +
  stat_slab(side = "right", show.legend = F,
            scale = 0.6, # defines the height that a slab can reach
            position = position_dodge(width=.6), # distance between elements for dodging
            aes(fill_ramp = stat(level), fill=as.factor(F99o)), 
            .width = c(.50,.95,1)) +  # set shading
  stat_summary(fun.data = "mean_cl_normal", show.legend = F, size = .4,
               position = position_dodge2nudge(x=.05,width = .8)) +
# styling
  scale_fill_ramp_discrete(from='black', aesthetics = "fill_ramp")+ # set ramping color
  guides( # change name and display of legend elements
          color="none") + # suppresses color legend item) 
  scale_colour_manual(values = cbPalette, aesthetics = c("colour","fill"))+
  theme(plot.caption = element_text(face = "italic")) +
  xlab(str_wrap("Svarskategorier för item F99", width = 40)) +
  ylab("Indexvärde") +
  ggtitle("Jämförelse av ett item och indexvärde") +
  labs(caption = "Punkterna visar medelvärden för indexvärde, fälten motsvarar 50% och 95% av distributionen.",
       subtitle = "Item F99o - Hur ofta tycker du att det är riktigt härligt att leva?") +
  coord_flip()

4.22 Jämförelser befintligt och nytt index

Rapporter utifrån Stockholmsenkäten har hittills räknat samman ett index för “Psykisk hälsa”, baserat på enkätsvar för de fem items som i denna analys återfinns under rubriken “Psykosomatiska besvär”

Code

items.psm <- c("F88","F91","F93","F95","F98")

df.index <- df.all %>%
  select(items.psm, Kön, ARSKURS, ar, SkolSDO)

RIlistitems(df.index)

itemnr	item
F88	Hur ofta har du haft huvudvärk detta läsår?
F91	Hur ofta har du dålig aptit?
F93	Hur ofta har du under detta läsår haft ”nervös mage” (t.ex. magknip, magkramper, orolig mage, illamående, gaser, förstoppning eller diarré)?
F95	Hur ofta har du under detta läsår haft svårt att somna?
F98	Hur ofta har det hänt under detta läsår att du sovit oroligt och vaknat under natten?

Indexberäkningen som använts utgår från att svaren översätts till siffror och sedan beräknas ett medelvärde av alla items. Lägsta svarsalternativet (“Aldrig”) ersätts med siffran 100, näst lägsta med 75, och så vidare. Det innebär att ett högt värde motsvarar låg förekomst av besvär.

I Stockholm Stads rapport för år 2022 återges resultatet på detta vis:

Vi kan återskapa figuren utifrån data från 2020 för årskurs 9.

Code

df.indexR <- df.index %>%
  mutate(across(items.psm, ~ recode(.x, "0=100;1=75;2=50;3=25;4=0"))) %>%
  na.omit() %>%
  mutate(PSFindex = rowMeans(across(items.psm)))

library(ggchicklet)
# create a wide palette based on RISE colors
RISEpalette1 <- colorRampPalette(colors = c("#009ca6", "#e83c63", "#ffe500"))(6)
# "#009CA6" "#5C758B" "#B94F70" "#EC5D4F" "#F5A127" "#FFE500"
RISEpalette2 <- colorRampPalette(colors = c("#009ca6", "#e83c63", "#ffe500"))(8)

df.indexR %>%
  filter(ar == 2020) %>%
  filter(Kön %in% c("Pojke", "Flicka")) %>%
  filter(!is.na(SkolSDO)) %>%
  filter(ARSKURS == "Åk 9") %>% 
  group_by(SkolSDO, Kön) %>%
  summarise(PSFmean = round(mean(PSFindex, na.rm = T), 0)) %>%
  ggplot(aes(y = PSFmean, x = SkolSDO, fill = Kön, color = Kön)) +
  geom_chicklet(
    radius = grid::unit(1.5, "mm"),
    position = "dodge"
  ) +
  geom_hline(data = . %>% filter(Kön == "Pojke"), aes(yintercept = mean(PSFmean)), linetype = 2, color = "darkgrey") +
  geom_hline(data = . %>% filter(Kön == "Flicka"), aes(yintercept = mean(PSFmean)), linetype = 3, color = "black") +
  geom_text(aes(label = PSFmean, y = PSFmean + 2),
    position = position_dodge(width = .9),
    color = "black", size = 3.8
  ) +

  # annotate("text", label = "Medelvärde", y = 57, x = 3, color = "darkgrey") +
  scale_fill_manual("Kön", values = RISEpalette1[c(1, 5)], aesthetics = c("color", "fill")) +
  scale_x_discrete("Stadsdel", guide = guide_axis(n.dodge = 3)) +
  scale_y_continuous(limits = c(0,100)) +
  theme_minimal(base_size = 12) +
  theme_rise(axissize = 14, titlesize = 16) +
  theme(text = element_text(family = "Lato")) +
  ylab("Medelvärde") +
  labs(
    title = "Psykisk hälsa indexvärde",
    subtitle = "Befintligt index, 2020, åk 9",
    caption = "Streckade linjer betecknar medelvärde för samtliga stadsdelar"
  )

4.22.1 Nytt index (5 items)

Här utgår vi från samma fem items som i det befintliga indexet, men kodar om svarskategorierna som inte fungerade som tänkt, och estimerar fram ett indexvärde för varje individ utifrån Rasch-parametrarna. För jämförbarhet skalas sedan mätvärdet om till att omfånget 0-100.

Frågorna mäter i huvudsak besvär, vilket medfört att indexet uppkallats efter det. För jämförbarhet med det befintliga indexet “psykisk hälsa”, har mätvärden reverserats, så högre värden = lägre nivåer av besvär. Samma gäller indexet med 8 items.

Code

# read item parameters from previous analysis
itemParams.psm <- read_csv("PsyksomItemParams.csv") %>% 
  as.matrix()

# create id variable to re-join thetas
df.index <- df.index %>% 
  mutate(id = seq.int(nrow(.)))

# recode items to correct response categories
df.indexNew <- df.index %>% 
  select(all_of(c(items.psm,"id"))) %>% 
  na.omit()

df.indexNew$F91 <- recode(df.indexNew$F91, "2=1;3=2;4=2", as.factor = FALSE)
df.indexNew$F93 <- recode(df.indexNew$F93, "4=3", as.factor = FALSE)
df.indexNew$F98 <- recode(df.indexNew$F98, "4=3", as.factor = FALSE)
df.indexNew$F95 <- recode(df.indexNew$F95, "4=3;3=2;2=1;1=0", as.factor = FALSE)

#RItileplot(df.indexNew %>% select(!id))
#itemParams.psm
df.indexNew$PSMscore <- df.indexNew %>% 
  select(!id) %>% 
  RIestThetas2(cpu = 8, itemParams = itemParams.psm)
df.index <- df.indexNew %>% 
  select(id,PSMscore) %>% 
  left_join(df.index,df.indexNew, by = "id")

# reverse score to match old index
df.index$PSMscore <- df.index$PSMscore*-1

Code

# rescale to match 0-100 range
df.index$PSMscore100 <- round(scales::rescale(df.index$PSMscore, to = c(0,100)),0)

4.22.2 Nytt index (8 items)

Det index som både representerar psykosomatiska och psykiska besvär innehåller 8 frågor, och beräknas utifrån items efter korrigering av svarskategorier.

Code

newIndex <- read_csv("PSFitemnr.csv")

items.psnew <- newIndex$itemnr

df.index8 <- df.all %>%
  select(items.psnew, Kön, ARSKURS, ar, SkolSDO)

RIlistitems(df.index8)

itemnr	item
F88	Hur ofta har du haft huvudvärk detta läsår?
F91	Hur ofta har du dålig aptit?
F92	Hur mycket skulle du vilja ändra på dig själv?
F93	Hur ofta har du under detta läsår haft ”nervös mage” (t.ex. magknip, magkramper, orolig mage, illamående, gaser, förstoppning eller diarré)?
F95	Hur ofta har du under detta läsår haft svårt att somna?
F97	Känner du dig slö och olustig?
F98	Hur ofta har det hänt under detta läsår att du sovit oroligt och vaknat under natten?
F99	Hur ofta tycker du att det är riktigt härligt att leva?

Code

# read item parameters from previous analysis
itemParams.psnew <- read_csv("itemParameters.csv") %>% 
  as.matrix()

# create id variable to re-join thetas
df.index8 <- df.index8 %>% 
  mutate(id = seq.int(nrow(.)))

# recode items to correct response categories
df.indexNew8 <- df.index8 %>% 
  select(all_of(c(items.psnew,"id"))) %>% 
  na.omit()

rcat1 <- c("F91")
rcat2 <- c("F95", "F97")
rcat3 <- c("F92", "F93", "F98", "F99")

for (i in rcat1) {
  df.indexNew8[[i]] <- recode(df.indexNew8[[i]], "2=1;3=2;4=2", as.factor = FALSE)
}
for (i in rcat2) {
  df.indexNew8[[i]] <- recode(df.indexNew8[[i]], "1=0;2=1;3=2;4=3", as.factor = FALSE)
}
for (i in rcat3) {
  df.indexNew8[[i]] <- recode(df.indexNew8[[i]], "4=3", as.factor = FALSE)
}

#RItileplot(df.indexNew %>% select(!id))
#itemParams.psm
df.indexNew8$PS8score <- df.indexNew8 %>% 
  select(!id) %>% 
  RIestThetas2(cpu = 8, itemParams = itemParams.psnew)
df.index8 <- df.indexNew8 %>% 
  select(id,PS8score) %>% 
  left_join(df.index8,df.indexNew8, by = "id")

# reverse score to match old index
df.index8$PS8score <- df.index8$PS8score*-1
# and rescaled version
df.index8$PS8score100 <- round(scales::rescale(df.index8$PS8score, to = c(0,100)),0)

Code

df.index %>% 
  filter(ar == 2020) %>% 
  filter(Kön %in% c("Pojke","Flicka")) %>% 
  filter(ARSKURS == "Åk 9") %>% 
  filter(!is.na(SkolSDO)) %>% 
  group_by(SkolSDO,Kön) %>% 
  summarise(PSFmean = round(mean(PSMscore100, na.rm = T),0)) %>% 
  ggplot(aes(y = PSFmean, x = SkolSDO, fill = Kön, color = Kön)) +
  geom_chicklet(
    radius = grid::unit(1.5, "mm"),
    position = "dodge"
  ) +
  geom_hline(data = . %>% filter(Kön == "Pojke"), aes(yintercept = mean(PSFmean)), linetype = 2, color = "darkgrey") +
  geom_hline(data = . %>% filter(Kön == "Flicka"), aes(yintercept = mean(PSFmean)), linetype = 3, color = "black") +
  geom_text(aes(label = PSFmean, y = PSFmean + 2),
    position = position_dodge(width = .9),
    color = "black", size = 3.8
  ) +

  # annotate("text", label = "Medelvärde", y = 57, x = 3, color = "darkgrey") +
  scale_fill_manual("Kön", values = RISEpalette1[c(1, 5)], aesthetics = c("color", "fill")) +
  scale_x_discrete("Stadsdel", guide = guide_axis(n.dodge = 3)) +
  scale_y_continuous(limits = c(0,100)) +
  theme_minimal(base_size = 12) +
  theme_rise(axissize = 14, titlesize = 16) +
  theme(text = element_text(family = "Lato")) +
  ylab("Medelvärde") +
  labs(title = "Psykosomatiska besvär indexvärde",
       subtitle = "5 items, år 2020, åk 9",
       caption = "Streckade linjer betecknar medelvärde per kön för samtliga stadsdelar")

Code

df.index8 %>% 
  filter(ar == 2020) %>% 
  filter(Kön %in% c("Pojke","Flicka")) %>% 
  filter(ARSKURS == "Åk 9") %>% 
  filter(!is.na(SkolSDO)) %>% 
  group_by(SkolSDO,Kön) %>% 
  summarise(PSFmean = round(mean(PS8score100, na.rm = T),0)) %>% 
  ggplot(aes(y = PSFmean, x = SkolSDO, fill = Kön, color = Kön)) +
  geom_chicklet(
    radius = grid::unit(1.5, "mm"),
    position = "dodge"
  ) +
  geom_hline(data = . %>% filter(Kön == "Pojke"), aes(yintercept = mean(PSFmean)), linetype = 2, color = "darkgrey") +
  geom_hline(data = . %>% filter(Kön == "Flicka"), aes(yintercept = mean(PSFmean)), linetype = 3, color = "black") +
  geom_text(aes(label = PSFmean, y = PSFmean + 2),
    position = position_dodge(width = .9),
    color = "black", size = 3.8
  ) +

  # annotate("text", label = "Medelvärde", y = 57, x = 3, color = "darkgrey") +
  scale_fill_manual("Kön", values = RISEpalette1[c(1, 5)], aesthetics = c("color", "fill")) +
  scale_x_discrete("Stadsdel", guide = guide_axis(n.dodge = 3)) +
  scale_y_continuous(limits = c(0,100)) +
  theme_minimal(base_size = 12) +
  theme_rise(axissize = 14, titlesize = 16) +
  theme(text = element_text(family = "Lato")) +
  ylab("Medelvärde") +
  labs(title = "Psykiska/psykosomatiska besvär indexvärde",
       subtitle = "8 items, år 2020, åk 9",
       caption = "Streckade linjer betecknar medelvärde per kön för samtliga stadsdelar")

4.22.3 Indexvärde över tid, åk 9

Code

årtal <- c(2006,2008,2010,2012,2014,2016,2018,2020,2022)
# library(ggiraph)
# psfTid <- df.indexR %>%
#   filter(Kön %in% c("Pojke", "Flicka")) %>%
#   filter(!is.na(SkolSDO)) %>%
#   group_by(ar, SkolSDO, Kön) %>%
#   summarise(PSFmean = round(mean(PSFindex, na.rm = T), 0),
#             stdev = round(sd(PSFindex, na.rm = T))) %>%
#   ggplot(aes(y = PSFmean, x = ar, color = Kön, group = Kön)) +
#   geom_line(linewidth = 1) +
#   geom_point_interactive(aes(tooltip = PSFmean),
#                          size = 1.5) +
#   geom_ribbon(aes(ymin = PSFmean-stdev, ymax = PSFmean+stdev, fill = Kön), 
#               alpha = 0.15, linetype = 0) +
#   geom_hline(data = . %>% filter(Kön == "Pojke"), aes(yintercept = mean(PSFmean)), 
#              linetype = 2, color = "black", alpha = 0.8) +
#   geom_hline(data = . %>% filter(Kön == "Flicka"), aes(yintercept = mean(PSFmean)), 
#              linetype = 3, color = "black", alpha = 0.8) +
#   scale_fill_manual("Kön", values = RISEpalette1[c(1, 5)], aesthetics = c("color", "fill")) +
#   scale_x_continuous("År",guide = guide_axis(n.dodge = 2), breaks = årtal) +
#   scale_y_continuous("Medelvärde", limits = c(0,100)) +
#   labs(title = "Psykisk ohälsa indexvärde",
#        caption = "Skuggade fält indikerar en standardavvikelse") +
#   facet_wrap(~SkolSDO) +
#   theme_minimal(base_size = 11) +
#   theme_rise()
# 
# ggiraph(ggobj = psfTid)

df.indexR %>%
  filter(Kön %in% c("Pojke", "Flicka")) %>%
  filter(!is.na(SkolSDO)) %>%
  filter(ARSKURS == "Åk 9") %>% 
  group_by(ar, SkolSDO, Kön) %>%
  summarise(PSFmean = round(mean(PSFindex, na.rm = T), 0),
            stdev = round(sd(PSFindex, na.rm = T))) %>%
  ggplot(aes(y = PSFmean, x = ar, color = Kön, group = Kön)) +
  geom_line(linewidth = 1) +
  geom_point(size = 2.5) +
  geom_ribbon(aes(ymin = PSFmean-stdev, ymax = PSFmean+stdev, fill = Kön), 
              alpha = 0.15, linetype = 0) +
  geom_hline(data = . %>% filter(Kön == "Pojke"), aes(yintercept = mean(PSFmean)), 
             linetype = 2, color = "black", alpha = 0.8) +
  geom_hline(data = . %>% filter(Kön == "Flicka"), aes(yintercept = mean(PSFmean)), 
             linetype = 3, color = "black", alpha = 0.8) +
  scale_fill_manual("Kön", values = RISEpalette1[c(1, 5)], aesthetics = c("color", "fill")) +
  scale_x_continuous("År",guide = guide_axis(n.dodge = 2), breaks = årtal) +
  scale_y_continuous("Medelvärde", limits = c(0,100)) +
  labs(title = "Psykisk hälsa indexvärde",
        subtitle = "5 items, befintligt index",
       caption = "Skuggade fält indikerar en standardavvikelse") +
  facet_wrap(~SkolSDO) +
  theme_minimal(base_size = 11) +
  theme_rise() +
  theme(text = element_text(family = "Lato"))

Code

# psfTidscore <- df.index %>%
#   filter(Kön %in% c("Pojke", "Flicka")) %>%
#   filter(!is.na(SkolSDO)) %>%
#   group_by(ar, SkolSDO, Kön) %>%
#   summarise(PSFmean = round(mean(PSMscore100, na.rm = T), 0),
#             stdev = round(sd(PSMscore100, na.rm = T))) %>%
#   ggplot(aes(y = PSFmean, x = ar, color = Kön, group = Kön)) +
#   geom_line(linewidth = 1) +
#   geom_point_interactive(aes(tooltip = PSFmean),
#                          size = 1.5) +
#   geom_ribbon(aes(ymin = PSFmean-stdev, ymax = PSFmean+stdev, fill = Kön), 
#               alpha = 0.15, linetype = 0) +
#   geom_hline(data = . %>% filter(Kön == "Pojke"), aes(yintercept = mean(PSFmean)), 
#              linetype = 2, color = "black", alpha = 0.8) +
#   geom_hline(data = . %>% filter(Kön == "Flicka"), aes(yintercept = mean(PSFmean)), 
#              linetype = 3, color = "black", alpha = 0.8) +
#   scale_fill_manual("Kön", values = RISEpalette1[c(1, 5)], aesthetics = c("color", "fill")) +
#   scale_x_continuous("År",guide = guide_axis(n.dodge = 2), breaks = årtal) +
#   scale_y_continuous("Medelvärde", limits = c(0,100)) +
#   labs(title = "Psykosomatiska besvär indexvärde",
#        caption = "Skuggade fält indikerar en standardavvikelse") +
#   facet_wrap(~SkolSDO) +
#   theme_minimal(base_size = 11) +
#   theme_rise()
# 
# ggiraph(ggobj = psfTidscore)

df.index %>%
  filter(Kön %in% c("Pojke", "Flicka")) %>%
  filter(!is.na(SkolSDO)) %>%
  filter(ARSKURS == "Åk 9") %>% 
  group_by(ar, SkolSDO, Kön) %>%
  summarise(PSFmean = round(mean(PSMscore100, na.rm = T), 0),
            stdev = round(sd(PSMscore100, na.rm = T))) %>%
  ggplot(aes(y = PSFmean, x = ar, color = Kön, group = Kön)) +
  geom_line(linewidth = 1) +
  geom_point(size = 2.5) +
  geom_ribbon(aes(ymin = PSFmean-stdev, ymax = PSFmean+stdev, fill = Kön), 
              alpha = 0.15, linetype = 0) +
  geom_hline(data = . %>% filter(Kön == "Pojke"), aes(yintercept = mean(PSFmean)), 
             linetype = 2, color = "black", alpha = 0.8) +
  geom_hline(data = . %>% filter(Kön == "Flicka"), aes(yintercept = mean(PSFmean)), 
             linetype = 3, color = "black", alpha = 0.8) +
  scale_fill_manual("Kön", values = RISEpalette1[c(1, 5)], aesthetics = c("color", "fill")) +
  scale_x_continuous("År",guide = guide_axis(n.dodge = 2), breaks = årtal) +
  scale_y_continuous("Medelvärde", limits = c(0,100)) +
  labs(title = "Psykosomatiska besvär indexvärde",
              subtitle = "5 items, nytt index",
       caption = "Skuggade fält indikerar en standardavvikelse") +
  facet_wrap(~SkolSDO) +
  theme_minimal(base_size = 11) +
  theme_rise() +
  theme(text = element_text(family = "Lato"))

Code

df.index8 %>%
  filter(Kön %in% c("Pojke", "Flicka")) %>%
  filter(!is.na(SkolSDO)) %>%
  filter(ARSKURS == "Åk 9") %>% 
  group_by(ar, SkolSDO, Kön) %>%
  summarise(PSFmean = round(mean(PS8score100, na.rm = T), 0),
            stdev = round(sd(PS8score100, na.rm = T))) %>%
  ggplot(aes(y = PSFmean, x = ar, color = Kön, group = Kön)) +
  geom_line(linewidth = 1) +
  geom_point(size = 2.5) +
  geom_ribbon(aes(ymin = PSFmean-stdev, ymax = PSFmean+stdev, fill = Kön), 
              alpha = 0.15, linetype = 0) +
  geom_hline(data = . %>% filter(Kön == "Pojke"), aes(yintercept = mean(PSFmean)), 
             linetype = 2, color = "black", alpha = 0.8) +
  geom_hline(data = . %>% filter(Kön == "Flicka"), aes(yintercept = mean(PSFmean)), 
             linetype = 3, color = "black", alpha = 0.8) +
  scale_fill_manual("Kön", values = RISEpalette1[c(1, 5)], aesthetics = c("color", "fill")) +
  scale_x_continuous("År",guide = guide_axis(n.dodge = 2), breaks = årtal) +
  scale_y_continuous("Medelvärde", limits = c(0,100)) +
  labs(title = "Psykiska/psykosomatiska besvär indexvärde",
              subtitle = "8 items, nytt index",
       caption = "Skuggade fält indikerar en standardavvikelse") +
  facet_wrap(~SkolSDO) +
  theme_minimal(base_size = 11) +
  theme_rise() +
  theme(text = element_text(family = "Lato"))

4.22.4 Fördelningar

Code

df.indexR %>% 
  filter(ar == 2020) %>% 
  filter(ARSKURS == "Åk 9") %>% 
  select(PSFindex) %>% 
  summary()

    PSFindex     
 Min.   :  0.00  
 1st Qu.: 40.00  
 Median : 55.00  
 Mean   : 55.49  
 3rd Qu.: 75.00  
 Max.   :100.00  

Code

df.indexR %>% 
  filter(ar == 2020) %>% 
  filter(ARSKURS == "Åk 9") %>% 
  select(PSFindex) %>% 
  pull() %>% 
  hist(., col = RISEprimGreen, main = "Psykisk ohälsa indexvärden")

Code

summary(df.indexR$PSFindex)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   0.00   40.00   60.00   57.54   75.00  100.00 

Code

df.index %>% 
  filter(ar == 2020) %>% 
  filter(ARSKURS == "Åk 9") %>% 
  select(PSMscore100) %>% 
  summary()

  PSMscore100    
 Min.   :  0.00  
 1st Qu.: 41.00  
 Median : 51.00  
 Mean   : 50.29  
 3rd Qu.: 62.00  
 Max.   :100.00  

Code

df.index %>% 
  filter(ar == 2020) %>%  
  filter(ARSKURS == "Åk 9") %>% 
  select(PSMscore100) %>% 
  pull() %>% 
  hist(., col = RISEprimGreen, main = "Psykosomatiska besvär 5 item")

Code

df.index8 %>% 
  filter(ar == 2020) %>% 
  filter(ARSKURS == "Åk 9") %>% 
  select(PS8score100) %>% 
  summary()

  PS8score100    
 Min.   :  0.00  
 1st Qu.: 43.00  
 Median : 50.00  
 Mean   : 50.87  
 3rd Qu.: 58.00  
 Max.   :100.00  

Code

df.index8 %>% 
  filter(ar == 2020) %>%
  filter(ARSKURS == "Åk 9") %>% 
  select(PS8score100) %>% 
  pull() %>% 
  hist(., col = RISEprimGreen, main = "Psykosomatiska besvär 8 item")

4.22.5 Violin- & boxplot jämförelse

Code

library(PrettyCols)
df.indexR %>%
  filter(Kön %in% c("Pojke", "Flicka")) %>%
  filter(!is.na(SkolSDO)) %>%
  filter(ARSKURS == "Åk 9") %>%
  mutate(År = factor(ar)) %>% 
  ggplot(aes(x = År, y = PSFindex, group = År)) +
  geom_jitter(shape = 16, position = position_jitter(0.3), alpha = 0.05) +
  geom_violin(aes(fill = År), alpha= 0.85) +
  geom_boxplot(width = .3, outlier.shape = NA, alpha = 0.3, notch = TRUE) +
  xlab("År") +
  ylab("Indexvärde") +
  # scale_color_manual(
  #   values = RISEpalette1[c(1, 5)],
  #   aesthetics = c("fill", "color")
  # ) +
   scale_fill_pretty_d(name = "TealGreens") +
  labs(
    title = "Psykisk hälsa",
    subtitle = "Befintligt index",
    caption = "I boxen ryms 50% av respondenterna. Strecket indikerar medianvärdet"
  ) +
  theme_minimal() +
  theme_rise() +
  theme(
    text = element_text(family = "Lato")
  )

Code

df.index %>%
  filter(Kön %in% c("Pojke", "Flicka")) %>%
  filter(!is.na(SkolSDO)) %>%
  filter(ARSKURS == "Åk 9") %>%
  mutate(År = factor(ar)) %>% 
  ggplot(aes(x = År, y = PSMscore100, group = År)) +
  geom_jitter(shape = 16, position = position_jitter(0.3), alpha = 0.05) +
  geom_violin(aes(fill = År), alpha= 0.85) +
  geom_boxplot(width = .3, outlier.shape = NA, alpha = 0.3, notch = TRUE) +
  xlab("År") +
  ylab("Indexvärde") +
  # scale_color_manual(
  #   values = RISEpalette1[c(1, 5)],
  #   aesthetics = c("fill", "color")
  # ) +
   scale_fill_pretty_d(name = "TealGreens") +
  labs(
    title = "Psykosomatiska besvär",
    subtitle = "Nytt 5-item index",
    caption = "I boxen ryms 50% av respondenterna. Strecket indikerar medianvärdet"
  ) +
  theme_minimal() +
  theme_rise() +
  theme(
    text = element_text(family = "Lato")
  )

Code

df.index8 %>%
  filter(Kön %in% c("Pojke", "Flicka")) %>%
  filter(!is.na(SkolSDO)) %>%
  filter(ARSKURS == "Åk 9") %>%
  mutate(År = factor(ar)) %>% 
  ggplot(aes(x = År, y = PS8score100, group = År)) +
  geom_jitter(shape = 16, position = position_jitter(0.3), alpha = 0.05) +
  geom_violin(aes(fill = År), alpha= 0.85) +
  geom_boxplot(width = .3, outlier.shape = NA, alpha = 0.3, notch = TRUE) +
  xlab("År") +
  ylab("Indexvärde") +
  # scale_color_manual(
  #   values = RISEpalette1[c(1, 5)],
  #   aesthetics = c("fill", "color")
  # ) +
   scale_fill_pretty_d(name = "TealGreens") +
  labs(
    title = "Psykiska/psykosomatiska besvär",
    subtitle = "Nytt 8-item index",
    caption = "I boxen ryms 50% av respondenterna. Strecket indikerar medianvärdet"
  ) +
  theme_minimal() +
  theme_rise() +
  theme(
    text = element_text(family = "Lato")
  )

4.22.5.1 Fördelat på kön

Code

library(PrettyCols)
df.indexR %>%
  filter(Kön %in% c("Pojke", "Flicka")) %>%
  filter(!is.na(SkolSDO)) %>%
  filter(ARSKURS == "Åk 9") %>%
  mutate(År = factor(ar)) %>% 
  ggplot(aes(x = År, y = PSFindex, group = År)) +
  geom_jitter(shape = 16, position = position_jitter(0.3), alpha = 0.05) +
  geom_violin(aes(fill = År), alpha= 0.85) +
  geom_boxplot(width = .3, outlier.shape = NA, alpha = 0.3, notch = TRUE) +
  xlab("År") +
  ylab("Indexvärde") +
  # scale_color_manual(
  #   values = RISEpalette1[c(1, 5)],
  #   aesthetics = c("fill", "color")
  # ) +
   scale_fill_pretty_d(name = "Rainbow") +
  labs(
    title = "Psykisk hälsa",
    subtitle = "Befintligt index",
    caption = "I boxen ryms 50% av respondenterna. Strecket indikerar medianvärdet"
  ) +
  theme_minimal() +
  theme_rise() +
  theme(
    text = element_text(family = "Lato")
  ) + 
  facet_wrap(~Kön)

Code

df.index %>%
  filter(Kön %in% c("Pojke", "Flicka")) %>%
  filter(!is.na(SkolSDO)) %>%
  filter(ARSKURS == "Åk 9") %>%
  mutate(År = factor(ar)) %>% 
  ggplot(aes(x = År, y = PSMscore100, group = År)) +
  geom_jitter(shape = 16, position = position_jitter(0.3), alpha = 0.05) +
  geom_violin(aes(fill = År), alpha= 0.85) +
  geom_boxplot(width = .3, outlier.shape = NA, alpha = 0.3, notch = TRUE) +
  xlab("År") +
  ylab("Indexvärde") +
  # scale_color_manual(
  #   values = RISEpalette1[c(1, 5)],
  #   aesthetics = c("fill", "color")
  # ) +
   scale_fill_pretty_d(name = "Rainbow") +
  labs(
    title = "Psykosomatiska besvär",
    subtitle = "Nytt 5-item index",
    caption = "I boxen ryms 50% av respondenterna. Strecket indikerar medianvärdet"
  ) +
  theme_minimal() +
  theme_rise() +
  theme(
    text = element_text(family = "Lato")
  ) +
  facet_wrap(~Kön)

Code

df.index8 %>%
  filter(Kön %in% c("Pojke", "Flicka")) %>%
  filter(!is.na(SkolSDO)) %>%
  filter(ARSKURS == "Åk 9") %>%
  mutate(År = factor(ar)) %>% 
  ggplot(aes(x = År, y = PS8score100, group = År)) +
  geom_jitter(shape = 16, position = position_jitter(0.3), alpha = 0.05) +
  geom_violin(aes(fill = År), alpha= 0.85) +
  geom_boxplot(width = .3, outlier.shape = NA, alpha = 0.3, notch = TRUE) +
  xlab("År") +
  ylab("Indexvärde") +
  # scale_color_manual(
  #   values = RISEpalette1[c(1, 5)],
  #   aesthetics = c("fill", "color")
  # ) +
   scale_fill_pretty_d(name = "Rainbow") +
  labs(
    title = "Psykiska/psykosomatiska besvär",
    subtitle = "Nytt 8-item index",
    caption = "I boxen ryms 50% av respondenterna. Strecket indikerar medianvärdet"
  ) +
  theme_minimal() +
  theme_rise() +
  theme(
    text = element_text(family = "Lato")
  ) + 
  facet_wrap(~Kön)

4.22.6 Tak-/golveffekter

Code

df.indexR %>% 
  group_by(PSFindex) %>% 
  reframe(n = n()) %>% 
  mutate(Procent = round(100*n/sum(n),1)) %>% 
  formattable(., table.attr = 'class=\"table table-striped\" style="font-size: 15px; font-family: Lato; width: 50%"')

PSFindex	n	Procent
0	577	0.7
5	865	1.0
10	1213	1.4
15	1654	1.9
20	2332	2.6
25	2937	3.3
30	3580	4.0
35	4307	4.9
40	5332	6.0
45	5841	6.6
50	6639	7.5
55	7077	8.0
60	7320	8.3
65	7190	8.1
70	7028	7.9
75	6367	7.2
80	5752	6.5
85	4339	4.9
90	3362	3.8
95	2271	2.6
100	2537	2.9

Code

df.index %>% 
  group_by(PSMscore100) %>% 
  reframe(n = n()) %>% 
  mutate(Procent = round(100*n/sum(n),1)) %>% 
  formattable(., table.attr = 'class=\"table table-striped\" style="font-size: 15px; font-family: Lato; width: 50%"')

PSMscore100	n	Procent
0	1958	2.0
18	2503	2.5
26	3512	3.6
32	4807	4.9
37	5733	5.8
41	6804	6.9
44	8239	8.4
48	9286	9.4
51	9385	9.5
54	9501	9.6
58	9154	9.3
62	8324	8.5
66	7015	7.1
72	5374	5.5
81	3564	3.6
100	3307	3.4

Code

df.index8 %>% 
  group_by(PS8score100) %>% 
  reframe(n = n()) %>% 
  mutate(Procent = round(100*n/sum(n),1)) %>% 
  formattable(., table.attr = 'class=\"table table-striped\" style="font-size: 15px; font-family: Lato; width: 50%"')

PS8score100	n	Procent
0	338	0.4
15	572	0.6
22	876	0.9
27	1275	1.3
31	1753	1.8
34	2293	2.4
36	2832	3.0
39	3543	3.7
41	3954	4.1
43	4712	4.9
45	5401	5.6
47	5822	6.1
48	6385	6.7
50	6624	6.9
52	6834	7.1
54	6938	7.2
56	6521	6.8
58	6218	6.5
61	5801	6.0
64	5021	5.2
67	4096	4.3
71	3373	3.5
76	2297	2.4
84	1467	1.5
100	1008	1.1

4.23 Programvara som använts för analyserna

Code

pkgs <- cite_packages(cite.tidyverse = TRUE, 
                      output = "table",
                      bib.file = "grateful-refs.bib",
                      include.RStudio = TRUE)
formattable(pkgs, 
            table.attr = 'class=\"table table-striped\" style="font-size: 15px; font-family: Lato; width: 80%"')

Package	Version	Citation
arrow	10.0.0	Richardson et al. (2022)
base	4.2.2	R Core Team (2022)
car	3.1.1	Fox and Weisberg (2019)
catR	3.17	Magis and Raîche (2012); Magis and Barrada (2017)
colorspace	2.1.0	Zeileis, Hornik, and Murrell (2009); Stauffer et al. (2009); Zeileis et al. (2020)
cowplot	1.1.1	Wilke (2020)
eRm	1.0.2	Mair and Hatzinger (2007b); Mair and Hatzinger (2007a); Hatzinger and Rusch (2009); Rusch, Maier, and Hatzinger (2013); Koller, Maier, and Hatzinger (2015); Debelak and Koller (2019); Mair, Hatzinger, and Maier (2021)
extrafont	0.18	Chang (2022)
foreach	1.5.2	Microsoft and Weston (2022)
formattable	0.2.1	Ren and Russell (2021)
furrr	0.3.1	Vaughan and Dancho (2022)
ggchicklet	0.6.0	Rudis (2022)
ggdist	3.2.0	Kay (2022)
gghalves	0.1.4	Tiedemann (2022)
ggpp	0.4.5	Aphalo (2022)
ggrepel	0.9.3	Slowikowski (2023)
glue	1.6.2	Hester and Bryan (2022)
grateful	0.1.11	Rodríguez-Sánchez, Jackson, and Hutchins (2022)
kableExtra	1.3.4	Zhu (2021)
knitr	1.42	Xie (2014); Xie (2015); Xie (2023)
matrixStats	0.63.0	Bengtsson (2022)
mirt	1.37.1	Chalmers (2012)
PrettyCols	1.0.1	Rennie (2023)
psych	2.2.9	Revelle (2022)
psychotree	0.16.0	Trepte and Verbeet (2010); Strobl, Wickelmaier, and Zeileis (2011); Strobl, Kopf, and Zeileis (2015); Komboz, Zeileis, and Strobl (2018); Wickelmaier and Zeileis (2018)
reshape	0.8.9	Wickham (2007)
RISEkbmRasch	0.1.11.1	Johansson (2023)
rmarkdown	2.20	Xie, Allaire, and Grolemund (2018); Xie, Dervieux, and Riederer (2020); Allaire et al. (2023)
scales	1.2.1	Wickham and Seidel (2022)
tidyverse	2.0.0	Wickham et al. (2019)

4.24 Referenser

Allaire, JJ, Yihui Xie, Jonathan McPherson, Javier Luraschi, Kevin Ushey, Aron Atkins, Hadley Wickham, Joe Cheng, Winston Chang, and Richard Iannone. 2023. Rmarkdown: Dynamic Documents for r. https://github.com/rstudio/rmarkdown.

Aphalo, Pedro J. 2022. Ggpp: Grammar Extensions to ’Ggplot2’. https://CRAN.R-project.org/package=ggpp.

Bengtsson, Henrik. 2022. matrixStats: Functions That Apply to Rows and Columns of Matrices (and to Vectors). https://CRAN.R-project.org/package=matrixStats.

Chalmers, R. Philip. 2012. “mirt: A Multidimensional Item Response Theory Package for the R Environment.” Journal of Statistical Software 48 (6): 1–29. https://doi.org/10.18637/jss.v048.i06.

Chang, Winston. 2022. Extrafont: Tools for Using Fonts. https://CRAN.R-project.org/package=extrafont.

Debelak, Rudolf, and Ingrid Koller. 2019. “Testing the Local Independence Assumption of the Rasch Model With Q3-Based Nonparametric Model Tests.” Applied Psychological Measurement. https://doi.org/10.1177/0146621619835501.

Fox, John, and Sanford Weisberg. 2019. An R Companion to Applied Regression. Third. Thousand Oaks CA: Sage. https://socialsciences.mcmaster.ca/jfox/Books/Companion/.

Hatzinger, Reinhold, and Thomas Rusch. 2009. “IRT models with relaxed assumptions in eRm: A manual-like instruction.” Psychology Science Quarterly 51.

Hester, Jim, and Jennifer Bryan. 2022. Glue: Interpreted String Literals. https://CRAN.R-project.org/package=glue.

Johansson, Magnus. 2023. RISEkbmRasch: Psychometric Analysis in r with Rasch Measurement Theory. https://github.com/pgmj/RISEkbmRasch.

Kay, Matthew. 2022. ggdist: Visualizations of Distributions and Uncertainty. https://doi.org/10.5281/zenodo.3879620.

Koller, Ingrid, Marco Johannes Maier, and Reinhold Hatzinger. 2015. “An Empirical Power Analysis of Quasi-Exact Tests for the Rasch Model: Measurement Invariance in Small Samples.” Methodology 11. https://doi.org/10.1027/1614-2241/a000090.

Komboz, Basil, Achim Zeileis, and Carolin Strobl. 2018. “Tree-Based Global Model Tests for Polytomous Rasch Models.” Educational and Psychological Measurement 78 (1): 128–66. https://doi.org/10.1177/0013164416664394.

Magis, David, and Juan Ramon Barrada. 2017. “Computerized Adaptive Testing with R: Recent Updates of the Package catR.” Journal of Statistical Software, Code Snippets 76 (1): 1–19. https://doi.org/10.18637/jss.v076.c01.

Magis, David, and Gilles Raîche. 2012. “Random Generation of Response Patterns Under Computerized Adaptive Testing with the R Package catR.” Journal of Statistical Software 48 (8): 1–31. https://doi.org/10.18637/jss.v048.i08.

Mair, Patrick, and Reinhold Hatzinger. 2007a. “CML based estimation of extended Rasch models with the eRm package in R.” Psychology Science 49.

———. 2007b. “Extended Rasch modeling: The eRm package for the application of IRT models in R.” Journal of Statistical Software 20. https://www.jstatsoft.org/v20/i09.

Mair, Patrick, Reinhold Hatzinger, and Marco Johannes Maier. 2021. eRm: Extended Rasch Modeling. https://cran.r-project.org/package=eRm.

Microsoft, and Steve Weston. 2022. Foreach: Provides Foreach Looping Construct. https://CRAN.R-project.org/package=foreach.

R Core Team. 2022. R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing. https://www.R-project.org/.

Ren, Kun, and Kenton Russell. 2021. Formattable: Create ’Formattable’ Data Structures. https://CRAN.R-project.org/package=formattable.

Rennie, Nicola. 2023. PrettyCols: Pretty Colour Palettes. https://CRAN.R-project.org/package=PrettyCols.

Revelle, William. 2022. Psych: Procedures for Psychological, Psychometric, and Personality Research. Evanston, Illinois: Northwestern University. https://CRAN.R-project.org/package=psych.

Richardson, Neal, Ian Cook, Nic Crane, Dewey Dunnington, Romain François, Jonathan Keane, Dragoș Moldovan-Grünfeld, Jeroen Ooms, and Apache Arrow. 2022. Arrow: Integration to ’Apache’ ’Arrow’. https://CRAN.R-project.org/package=arrow.

Rodríguez-Sánchez, Francisco, Connor P. Jackson, and Shaurita D. Hutchins. 2022. Grateful: Facilitate Citation of r Packages. https://github.com/Pakillo/grateful.

Rudis, Bob. 2022. Ggchicklet: Create ’Chicklet’ (Rounded Segmented Column) Charts. https://git.rud.is/hrbrmstr/ggchicklet.

Rusch, Thomas, Marco Johannes Maier, and Reinhold Hatzinger. 2013. “Linear logistic models with relaxed assumptions in R.” In Algorithms from and for Nature and Life, edited by Berthold Lausen, Dirk van den Poel, and Alfred Ultsch. Studies in Classification, Data Analysis, and Knowledge Organization. New York: Springer. https://doi.org/10.1007/978-3-319-00035-0_34.

Slowikowski, Kamil. 2023. Ggrepel: Automatically Position Non-Overlapping Text Labels with ’Ggplot2’. https://CRAN.R-project.org/package=ggrepel.

Stauffer, Reto, Georg J. Mayr, Markus Dabernig, and Achim Zeileis. 2009. “Somewhere over the Rainbow: How to Make Effective Use of Colors in Meteorological Visualizations.” Bulletin of the American Meteorological Society 96 (2): 203–16. https://doi.org/10.1175/BAMS-D-13-00155.1.

Strobl, Carolin, Julia Kopf, and Achim Zeileis. 2015. “Rasch Trees: A New Method for Detecting Differential Item Functioning in the Rasch Model.” Psychometrika 80 (2): 289–316. https://doi.org/10.1007/s11336-013-9388-3.

Strobl, Carolin, Florian Wickelmaier, and Achim Zeileis. 2011. “Accounting for Individual Differences in Bradley-Terry Models by Means of Recursive Partitioning.” Journal of Educational and Behavioral Statistics 36 (2): 135–53. https://doi.org/10.3102/1076998609359791.

Tiedemann, Frederik. 2022. Gghalves: Compose Half-Half Plots Using Your Favourite Geoms. https://CRAN.R-project.org/package=gghalves.

Trepte, Sabine, and Markus Verbeet, eds. 2010. Allgemeinbildung in Deutschland – Erkenntnisse Aus Dem SPIEGEL Studentenpisa-Test. Wiesbaden: VS Verlag.

Vaughan, Davis, and Matt Dancho. 2022. Furrr: Apply Mapping Functions in Parallel Using Futures. https://CRAN.R-project.org/package=furrr.

Wickelmaier, Florian, and Achim Zeileis. 2018. “Using Recursive Partitioning to Account for Parameter Heterogeneity in Multinomial Processing Tree Models.” Behavior Research Methods 50 (3): 1217–33. https://doi.org/10.3758/s13428-017-0937-z.

Wickham, Hadley. 2007. “Reshaping Data with the Reshape Package.” Journal of Statistical Software 21 (12). https://www.jstatsoft.org/v21/i12/.

Wickham, Hadley, Mara Averick, Jennifer Bryan, Winston Chang, Lucy D’Agostino McGowan, Romain François, Garrett Grolemund, et al. 2019. “Welcome to the tidyverse.” Journal of Open Source Software 4 (43): 1686. https://doi.org/10.21105/joss.01686.

Wickham, Hadley, and Dana Seidel. 2022. Scales: Scale Functions for Visualization. https://CRAN.R-project.org/package=scales.

Wilke, Claus O. 2020. Cowplot: Streamlined Plot Theme and Plot Annotations for ’Ggplot2’. https://CRAN.R-project.org/package=cowplot.

Xie, Yihui. 2014. “Knitr: A Comprehensive Tool for Reproducible Research in R.” In Implementing Reproducible Computational Research, edited by Victoria Stodden, Friedrich Leisch, and Roger D. Peng. Chapman; Hall/CRC.

———. 2015. Dynamic Documents with R and Knitr. 2nd ed. Boca Raton, Florida: Chapman; Hall/CRC. https://yihui.org/knitr/.

———. 2023. Knitr: A General-Purpose Package for Dynamic Report Generation in r. https://yihui.org/knitr/.

Xie, Yihui, J. J. Allaire, and Garrett Grolemund. 2018. R Markdown: The Definitive Guide. Boca Raton, Florida: Chapman; Hall/CRC. https://bookdown.org/yihui/rmarkdown.

Xie, Yihui, Christophe Dervieux, and Emily Riederer. 2020. R Markdown Cookbook. Boca Raton, Florida: Chapman; Hall/CRC. https://bookdown.org/yihui/rmarkdown-cookbook.

Zeileis, Achim, Jason C. Fisher, Kurt Hornik, Ross Ihaka, Claire D. McWhite, Paul Murrell, Reto Stauffer, and Claus O. Wilke. 2020. “colorspace: A Toolbox for Manipulating and Assessing Colors and Palettes.” Journal of Statistical Software 96 (1): 1–49. https://doi.org/10.18637/jss.v096.i01.

Zeileis, Achim, Kurt Hornik, and Paul Murrell. 2009. “Escaping RGBland: Selecting Colors for Statistical Graphics.” Computational Statistics & Data Analysis 53 (9): 3259–70. https://doi.org/10.1016/j.csda.2008.11.033.

Zhu, Hao. 2021. kableExtra: Construct Complex Table with ’Kable’ and Pipe Syntax. https://CRAN.R-project.org/package=kableExtra.